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>V Abstract 

We numerically study the static structure and the mechanical response of two-dimensional gran- 
ular piles. The piles consist of polydisperse disks with and without friction. Special attention is 
paid for the rigid grain limit by examining the systems with various disk elasticities. It is shown 
that the static pile structure of frictionless disks becomes isostatic in the rigid limit, while the 
isostaticity of frictional pile depends on the pile forming procedure, but in the case of the infinite 
friction is effective, the structure becomes very close to isostatic in the rigid limit. The mechanical 
response of the piles are studied by infinitesimally displacing one of the disks at the bottom. It is 
shown that the total amount of the displacement in the pile caused by the perturbation diverges 
in the case of frictionless pile as it becomes isostatic, while the response remains finite for the 
frictional pile. In the frictionless isostatic pile, the displacement response in each sample behaves 
■ rather complicated way, but its average shows wave like propagation. 
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I. INTRODUCTION 



A pile of granular material is often modeled as an assembly of rigid particles. This simple 
picture, however, causes some conceptual difficulties when one consider force distribution, 
or stress, in the pile. 

In an ordinary solid that consists of atoms, which are deformable elements, the stress 
inside is determined by the stress balance equations with the constitutive relation between 
the stress and strain. No deformation is allowed, however, in rigid elements, in which case 
the forces acting on each element cannot be determined uniquely for a given configuration 
in general. 

There is a special class of stable pile structures called isostatic, or marginally rigid: A 
pile structure is isostatic when the forces acting between the elements are uniquely deter- 
mined only from externally applied forces and the pile structure without any information 
on deformation of the elements. For such a pile, the total number of balance equations for 
the forces and torques acting each particle should be equal to the number of independent 
components of forces, which means the average coordination number z is 2d for frictionless 
spheres and d + 1 for frictional spheres or aspherical particles in <i-dimensions if the pile is 
isostatic. 

In a real pile the isostaticity is not always satisfied; in the overconstrained case, where 
the number of conditions is larger than the number of forces, the pile is unstable because 
there are no sets of force that satisfy the balance equations. On the other hand, in the 
underconstrained case, where the number of conditions is smaller than the number of forces, 
the forces cannot be determined uniquely from the macroscopic structure information of the 
pile; the macroscopic friction force depends on the piling history. It has been conjectured 
that a stable pile of frictionless rigid particles forms an isostatic structure Q]. 

The isostaticity of a pile structure has been tested numerically and experimentally for 
both the frictionless and frictional cases by counting the coordination numbers. Makse et 
a/.Q have performed numerical simulations for three dimensional sphere systems without 
the gravity and made compact aggregates of balls by compressing the systems by pushing 
the surrounding walls. By examining the coordination numbers in the zero-pressure limit, 
they concluded both the frictionless and frictional sphere systems become isostatic in the 
rigid limit. On the other hand, Silbert et a/.|2j have also performed numerical simulations on 
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the three-dimensional sphere system, but the way they made piles is different from that of 
Makse et al. They released the particles in the system at once under the gravity and waited 
until all the particles stop. Their conclusion is that the frictionless pile becomes isostatic 
in the rigid limit, but the structure of frictionalpile depends on the piling procedure and 
never becomes isostatic. Blumenfeld and Bali0, 0] have done simple table top experiments 
on the pile of two-dimensional non-circular grains made of cardboard. Piles are formed by 
collecting the grains scattered initially on a horizontal surface by sliding an open rectangular 
frame. They found that the higher the starting density of grains is, the more cooperative 
reconfiguration is taken place before they are stuck with each other, and the lower the ending 
pile density. They have concluded that the isostatic structure of frictional grain is achieved 
in the limiting case where the starting and ending density coincide. 

The isostatic structure, if it is realized in a real pile, should be reflected in the mechanical 
properties of the rigid granular pile. The mechanical properties of frictionless isostatic 
structure has been studied in some detail, and it has been shown that (i) the force chain 
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the force-force response and the displacement-displacement response pj], and (iii) the piles 



These features may correspond with some 



are very sensitive to an external perturbation 
properties of frictionless rigid granular pile. 

Regarding the mechanical response of the isostatic pile, using a simple lattice model of iso- 
static structure, MoukarzeljjJ has shown that the total response to an external perturbation 
diverges as the system approaches to the isostaticity due to the "pantograph effect" . As for 
the effects of friction, Moukarzel et al. js| have performed experiments on the two-dimensional 
rigid disk system and demonstrated that the displacement response function has the single 
peak Gaussian shape with diffusive broadening. They also performed computer simulations 
on the two-dimensional frictionless disk system and showed that the displacement response 
function have a double peaked shape with the wave like propagation, suggesting the friction 
plays an important role in the mechanical response. 

n 

In this paper, we present the detailed results of our numerical simulations 9] on the two- 
dimensional frictionless and frictional disks to show how the isostatic piles are formed for 
frictionless and frictional disks. We also investigate the mechanical response of the isostatic 
piles to an external perturbation, and demonstrate there exist important differences in the 
response between the frictionless and frictional piles when they approach the isostaticity. 
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After introducing the model and piling procedures in Sec. 2, we present the results for the 
static pile structure for frictionless and frictional disks in Sec. 3. The mechanical responses 
are investigated in Sec. 4. The summary of our results is given in Sec. 5. 



II. MODEL 

We perform molecular dynamics simulations on the system that consists of two- 
dimensional disks with linear elasticity and damping, which is usually called DEM(Discrete 
Element Method) in the engineer community. The system is polydisperse with a uniform 
distribution in the disk diameter over the range between 0.9<To and ctq, with <7o being the 
maximum diameter. The masses of the disks are assumed to be proportional to their areas: 
the mass of the disk with the diameter o"o is denoted by m®. The bottom of the system is 
made rough by attaching the disks with the interval <To, and we employ the periodic bound- 
ary condition in the horizontal direction. The number of disks N in the system is typically 
400, and the horizontal length of the system is 20<7o, thus we have approximately 20 layers 
of disks on average. 

Piles are formed by letting the system run under gravity with the acceleration g from 
initial configurations until all the disks stop moving. 

The disk at the position x(t) with mass m follows the Newtonian equation 

mx(t) = F(t), (1) 

where the force F(t) consists of the gravitational and contact forces from the neighboring 
disks in contact. We also use the "viscous equation" 

7i(f) = F(t) (2) 

for some cases to examine the effect of friction because the system that follows eq. (J2J) is 
stuck as soon as the force balance is achieved and is more affected by the friction. 

The two particles i and j at CEj(t) and Xj(t) with radii r$ and rj, respectively, are in 
contact when the overlap Sij given by 

§ij — Tj -\- Tj |ajj Xj\ (3) 

is positive. Then the particle j exerts the force Fij on the particle i 

F l3 =F% + F%, (4) 
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where F^ and F^ are the normal and tangential components of the force: 

F ij = Khjfiij -Wij, (5) 
F% = -hAsJij-jtvlj. (6) 

Here, hy and are the normal and tangential unit vectors, respectively; Asij is the tangen- 
tial displacement of the contact points after the contact, i>™ { v \j) is the normal (tangential) 
relative velocity, k n (k t ) is the normal (tangential) elastic constant, and 7„ (7^ is the normal 
(tangential) damping constant. Note that we assume no threshold for the disks to slip during 
the contact in eq.(jHJ), which corresponds to the case with the infinite friction coefficient. In 
the case of the frictionless disk, we simply set F\j = 0. 

In the actual simulations, we use 7„ = 2\fk^ [^/mo] for the frictionless case, and k t = 0.2k n 
and 7n = 7* = 2y/k^ [V m o] for the frictional case. In the simulations with the viscous 
equation (0), we take 7 = 5 [m a/^/^o] with 7 n = j t = in eqs.(jSJ) and ©• 

III. STATIC STRUCTURE OF PILES 

First, we study the isostaticity of granular pile formed through several procedures. 

Isostaticity is most easily checked by counting the average number of particles in contact 
with each particle, namely the average coordination number z. In order that the force acting 
between particle can be determined only by the contact network structure, the number 
of independent components of force should be equal to the number of force and torque 
balance equations. In the case of frictionless spherical grain, the number of independent 
force components is zN/2 because the contact forces have only radial component. The 
torque of each particle always balances and the number of force balance equation is Nd, 
therefore, we have z = 2d for the isostatic structure of frictionless spheres in ci-dimensions. 
In the case of frictional grain, we have z — d + 1 for both spherical and non-spherical grain, 
because the number for independent force components is zdN/2 and the number of the force 
and torque balance equations are dN and d(d — l)N/2, respectively. 

We perform molecular dynamics simulations to construct piles of frictionless and frictional 
disks using several procedures. We try two types of initial configurations: the triangular 
lattice and the random configuration (Fig. Q); the triangular lattice with the lattice constant 
ctq is not a regular lattice because the disks located at the lattice points are polydisperse. 
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The random configurations are prepared by randomly arranging disks with an area fraction 
of approximately 0.6. 

Simulations to form piles start from these configurations with zero particle velocity and 
finish when the kinetic energy of each disk becomes negligibly small, namely, smaller than 
lO~ 15 [m o g0"o]. The number of disks N is 400 and the system size is 20<7o in the horizontal 
direction, thus the number of layers in depth is 20 on average. We also try both eqs. 
and (J2J) for the time development. 

Figure Efa) shows the fc n -dependence of the coordination number z for the frictionless 
disks. It can be seen that the results does not depend on the preparation procedures very 
much and z converges to a number very close to 4 in the large k n limit for both initial 
configurations and time developments. The /c n -dependence is well represented by the power 
law 

z - z^ oc k~ a (7) 

as is shown in Fig. Of a). The parameters z^ and a are tabulated in Table |U 

As for the case of frictional disks, the results are shown in Fig. Efb). There are two 
things to be noted in comparison with the frictionless case: (i) the discrepancy among 
different preparation procedures is larger in the frictional case, and (ii) the limiting values 
of the coordination number are substantially different from 3, i.e. the value for the isostatic 
structure of the frictional grain in two dimensions(Table H}. 

There seems to be, however, the tendency that the limiting value of the coordination 
number z^ becomes closer to 3 in the case where the friction may produce more random pile 
configuration; namely from the Newtonian with the triangular lattice initial configurations 
to the viscous equation with the random initial configuration. 

These results should be compared with those by Silbert el al.{3\ They also constructed 
granular piles using DEM and concluded that the piles of frictionless spheres become isostatic 
in the rigid limit but that of frictional spheres does not become isostatic. 

The major difference between the present work and that of Silbert et al. lies in the follow- 
ing points: Silbert et a/.|3[ studied the three dimensional systems of mono-disperse spheres 
with a finite friction constant that follows the Newtonian equation while we investigate the 
two-dimensional system of polydisperse disks with the infinite friction constant that follows 
the Newtonian or the viscous equation. Both agree in the point that the pile structure in 
the frictional systems does not become isostatic in the same way as it does in the frictionless 
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system, but our results suggest that the isostaticity is achieved even for a frictional system 
in the certain limiting situation where the friction becomes very effective. 

To examine how this discrepancy arises, we plot the distribution of the ratio ( of the 
tangential force to the normal force at each contact. The comparison between the pile formed 
via the Newtonian equation and the pile via the viscous equation from the random initial 
configurations is given in Fig. H]for the disk elasticity k n = 10 6 [rriog / /<jq}, in which case 
the coordination numbers of the piles are z — 3.18 for the Newtonian equation and z = 3.06 
for the viscous equation. For the both cases we use the random initial configurations. One 
can see that the pile by the viscous equation contains more contacts with very large value of 
(, while the pile by the Newtonian equation has only contacts with ( smaller than 10 even 
though the infinite friction coefficient allows any value of (. 

If we use a finite value for the friction coefficient, some of the contacts in the pile via the 
viscous equation would slip to make the pile denser. This would results in a larger average 
coordination number, which means the pile structure becomes further away from isostatic. 

From these observations, we conclude that the infinite friction coefficient and the viscous 
equation for the time development in the pile forming process from a random initial config- 
uration makes the friction very effective, as a result, the pile becomes so decompacted that 
it becomes nearly isostatic. 

IV. MECHANICAL RESPONSE OF PILES 

Now, we study the mechanical response of the frictionless and frictional piles to the 
external perturbation and see how the response changes as the pile becomes closer to the 
isostatic. The perturbation is given by displacing one of the disks attached at the floor by 
<5r*o in the upward direction very slowly, and we observe the displacement of each disk 5ri 
caused in the pile by it. 

In order to examine properties of a given contact network, the size of the external dis- 
placement is taken small enough that the perturbation does not cause any change in the 
connectivity of the contact network in the pile. In the simulation, we take the external 
displacement as <5r*o = (0,mo/k n ), namely, the order of disk deformation. For such small 
external displacement, we have checked that the contact network in the pile does not change, 
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and Sri is proportional to 5r thus the relative displacement di defined by 



does not depend on \8tq\. 

The initial piles are prepared by the way described in the previous section with the 
random initial configuration and the Newtonian equation. The system size is 60<To in the 
horizontal direction and the number of disks N is 1200, thus the average number of layers 
in the depth is 20. 

Examples of pile response to the perturbation are shown in FigJSJ The filled disks in 
the bottom layer are fixed except one at the center marked with a red arrow, which disk 
is displaced upward. The directions and the distances of the displacements of the disks 
in the pile are shown by the colored arrows: the red, green and black arrows denotes the 
displacements \5di\ > 1, 0.5 < \5di\ < 1, and 0.1 < \5di\ < 0.5, respectively. The disks that 
move less than O.lSro are not marked. One can see the effects of the perturbation extend 
over long distance in the upper case, namely, the frictionless pile close to the isostaticity, 
while the effects decay within a short distance in the lower pile of frictional disks. 

A. Total Longitudinal Response of Displacement 

As a measure of response to the perturbation, we define the total response of displacement 
in the y-direction D y , or the total longitudinal response, as 



where 5yi is the y component of the displacement for the z-th disk. This quantity should be 
finite if the response is confined within a finite region, but can diverge in the infinite system 
if the response extend to the infinity. 

The results are shown in Figlfja), where the total response D y is plotted as a function of 
k n with which the examined piles are formed by the Newtonian equation from the random 
initial configuration. 

The marks 0's (H's) in FigEf a) denote the total response D y of the frictionless (frictional) 
disks in the frictionless (frictional) piles, respectively. 




(9) 
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In addition to these, we examine the frictional response to the perturbation in the friction- 
less piles (D's); namely, the response is calculated using the frictional interaction between 
disks, although the pile itself is prepared using the frictionless interaction, thus the pile 
structure has larger coordination numbers than that of real frictional piles with the same 
elastic constant k n . 

In FiglHfb), the same data are plotted against the coordination number z of the pile. 
One can see that the total response of the frictionless pile tends to diverge as z approaches 
4, while it remains finite for the frictional pile as z — > 3. It is interesting to see that the 
frictional responses of the frictionless piles are almost on the same curve with the frictional 
responses of the frictional piles when they are plotted against z, which suggests that the 
coordination number characterizes the mechanical response of the pile very well. 



B. Spatial Variation of Absolute Value of Longitudinal Response 

To understand the diverging total responses in the frictionless isostatic pile, we plot the 
averaged behavior of the spatial variation of absolute value of longitudinal response d y (r) 
defined by 

TV 

dy{r) =< \ d iy\ s ( r ~ r i) >, ( 10 ) 
1=1 

where < ... > denotes the statistical average. This quantity is related to D y by 




D y = / drd y {r). (11) 



In actual calculations, the spatial variation is calculated on a grid with the mesh spacing <Jo 
and the average is taken over a few hundreds realizations. 

The results are shown in FigJH where <i y (r)'s are plotted against x and y with the contour 
lines; The perturbation is applied at (x, y) = (0,0). 

We examine the two cases, namely, k n = 10 3 and 10 6 [m Q g/a }, for both the frictionless 
and frictional piles; the coordination number z = 4.86 and 3.98 for k n = 10 3 and 10 6 
[m g/(To\, respectively, for the frictionless pile, and z = 3.75 and 3.12 for the frictional pile; 
The piles with k n = 10 6 [mog/ao] are closer to the isostatic than those with k n = 10 3 
[mog/ao] for both the frictionless and frictional piles. 

Let us examine the frictionless case first. The response in the pile with k n = 10 3 [m g/a ], 
which is away from the isostatic, decays quickly as it departs from the point of perturbation. 



On the other hand, the situation is quite different for the pile with k n = 10 6 [mog/ao], which 
is close to the isostatic: the response does not decay along the y-axis with x — 0, and if 
one looks along the line parallel to the x-axis with constant y, one sees a plateau region 
where the response is constant. This plateau region, which seems to extend to the infinity, 
is responsible to the diverging behavior of the total response D y . 

On the other hand, in the case of frictional pile with k n = 10 6 [mog/ao], there is no 
tendency to develop the plateau region, although the response is larger than that for the 
pile of k n = 10 3 [mog/ao]. 

C. Averaged Displacement-Displacement Response Function 

Finally, we present the displacement-displacement response function A(r) for this exter- 
nal perturbation: 

N 

A(r)=<^d 4 «J(r-r i )>. (12) 

i=l 

The average is taken over a few hundreds realizations and the spatial dependence is calcu- 
lated on a grid with the mesh size ao- 

The displacement-displacement response function has been shown to be equal to the 
force-force response function for a frictionless isostatic structure j?]], but it should be noted 
that this correspondence does not hold in other cases. 

The results are shown in Figs. IHlandElfor the frictionless and frictional piles, respectively, 
for k n = 10 3 and 10 6 [mog/ao]. 

A x (r) is positive for x > and negative for x < while A y (r) is always positive when it 
is averaged. 

For both the frictionless and frictional piles, the region where A(r) ^ is larger for 
the pile whose coordination number z is smaller, namely, for the pile that is closer to the 
isostatic. 

It should be noted that the way A 3/ (r) extends is very different from that of d y (r), espe- 
cially in the frictionless pile close to the isostatic. The response function A y (r) propagates 
in the ^/-direction with the double peaked structure when one see it along the line parallel 
to the with constant y, while d y (r) develops the plateau region. 

This comparison shows that the response to the perturbation is not actually small in the 
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"low response region" of A y (r) between the two peaks for the frictionless isostatic pile, but 
the large response varies from sample to sample and they are averaged out to make A y (r) 
small. 

The way how this double peaked structure in A y (r) develops as the frictionless pile 
approaches to isostatic can be seen in the contour line plots in Fig^] (a). 

The response A(r) in the frictional pile shows smoother structure than that in the fric- 
tionless pile and does not seems to develop the double peaked structure as the pile becomes 
isostatic (Fig|HB(b)). 

V. SUMMARY 

We have investigated the structure and the mechanical response of the two-dimensional 
piles of disks with various disk elasticities; The piles are formed through a several deposition 
procedures under the gravity. 

As for the structure, we have shown the followings: The piles of frictionless disks become 
isostatic when the disks are very hard and they are not very sensitive to the preparation 
procedure, which is consistent with the conjecture that the pile of rigid grains is isostatic. 
On the other hand, for the piles of frictional disks with the infinite friction, the structure 
depends on the preparation process. If the pile is formed from a triangular lattice with 
the inertia following the Newtonian equation, the pile structure seems to be distinctively 
different from the isostatic one even in the rigid limit, as has been found in the previous 
work on the three-dimensional systemj^]. We have found, however, that the pile of frictional 
disks becomes very close to the isostatic one in the rigid grain limit when we employ the 
deposition process where the infinite friction is effective, namely, the viscous equation is used 
for the time development for the disks with the infinite friction constant from the random 
initial configurations. 

The role of friction for the frictional isostaticity is demonstrated by examining the dis- 
tribution of the ratio of the normal force to the tangential force; In the pile that is close to 
the isostaticity with frictional disks, the distribution of ( extends to very large value as of 
order 10 3 . This suggests that the frictional isostaticity is realized only in the cases where the 
exceptionally large friction coefficient is effective, and most of real stable piles with modest 
friction should be hyperstatic with history dependent forces even in the rigid limit. 
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We have also investigated the mechanical response of the piles of frictionless and frictional 
disks with a special attention on the isostaticity of the pile structure. We have examined 
the disk displacement caused by moving one of the disks at the bottom of the pile by an 
infinitesimally small distance. It has been shown that, for the frictionless pile with the 
isostatic structure, the response does not decay as it departs from the point of perturbation, 
and the total sum of the disk displacement diverges as the pile becomes isostatic, while the 
response decays quickly in the pile with structure far from isostaticity. It is found that the 
averaged longitudinal response A y becomes double peaked in the isostatic frictionless pile. 

As for the frictional pile, the response function always decays in a finite distance and the 
longitudinal response remains single peaked, which behavior does not change drastically as 
the structure approaches isostatic. 

The double peaked structure in the displacement-displacement response function for the 
frictionless isostatic pile may be compared with the hyperbolic stress propagation in the 
grannlar systemflQ due to the e q uivaW between the stress-stress response fonefon 
and the displacement-displacement response function in the frictionless isostatic pile|2J. The 
single peaked structure, on the other hand, may correspond to the diffusive stress propa- 
gation, but this correspondence is indirect because there is no equivalence between the 
displacement response and the stress response in the frictional or non-isostatic pile. 

It is interesting to note that, in the isostatic piles, the spatial distribution of averaged 
response of disk displacement is clearly different from that of absolute value of disk 
displacement; the former develops the double peaked structure while the latter shows the 
plateau structure, namely, in the region between the two peaks where the average displace- 
ment is small, the actual displacement response is not small in each sample, but just they 
are random and averaged out. This means that the averaged displacement response in the 
isostatic pile appears to propagate like a wave following a hyperbolic equation, but the way 
that displacement response propagates in each sample is quite complicated and does not 
look like a wave(Fig. EJ). 
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Eq. of motion 


Initial config. 


Frictionless 

Zoo a 


Frictional 

Zoo a 


Newtonian 


Triangular lattice 
Random config. 


3.97 0.68 

3.98 0.65 


3.15 0.49 
3.09 0.47 


Viscous 


Triangular lattice 
Random config. 


3.97 0.63 
3.97 0.64 


3.06 0.60 
3.04 0.46 



TABLE I: The limiting coordination numbers and the exponents a for various preparation 
procedures. 




FIG. 1: Two initial configurations: the triangular lattice with polydisperse disks (left) and the 
random configuration with a number density 0.6 (right). 
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(a) Frictionless Pile 
5.2 




Vkn 
(b) Frictional Pile 
4.4 



[oo/mog] 




10" 



10" J 10" z 
[ob/m g] 



FIG. 2: The coordination number z for various elastic constants k n for a frictionless pile (a) and 
frictional pile (b). The marks represent the pile preparation procedure: the Newtonian equation 
with the triangular lattice initial configuration (+), the Newtonian equation with the random initial 
configuration (x), the viscous equation with the triangular lattice initial configuration (*), and the 
viscous equation with the random initial configuration (□). Each mark represents the average of 
six to twelve realizations. 
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(a) Frictionless Pile 
1 




10 3 10 4 10 5 10 6 




[mog/oo] 



10 3 10 4 10 5 10 6 



[m g/o ] 



FIG. 3: The coordination numbers z — z^ for various elastic constants k n in the log- log scale for 
the frictionless pile (a) and the frictional pile (b). The same data are plotted using the same marks 
as those in Fig. [2] with Zqo listed in Table [I] 
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c 

FIG. 4: The distribution for the ratio ( of the tangential force to the normal force in the piles 
produced by the viscous equation (+) and that by Newtonian equation (□) from the random initial 
configurations with the disk elasticity k n = 10 6 [mo<7/co]- The average coordination numbers are 
z = 3.06 for the pile by the viscous equation and z = 3.18 for that by the Newtonian equation. 
Each plot represents average over about ten realizations. 
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FIG. 5: (Color online) The mechanical response caused by the small displacement of a disk at 
the bottom (denoted by an open (red) circle with a dark grey (red) arrow) in the frictionless pile 
(upper) and in the frictional pile (lower) with k n = 10 6 [mog/<7o]. The dark grey (red), light grey 
(green) , and black allows denote the displacement direction of the disks that move by the distance 
Sri > <5 r 0; 0.55rQ < 5ri < 5ro, and 0.1<5ro < 5ri < 0.5<5ro, respectively. 
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(b) 
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FIG. 6: The total longitudinal displacement D y caused by the perturbation in the frictionless pile 
(0), in the frictional pile (■), and for the frictional response in the piles formed through frictionless 
dynamics (□, see text), (a) D y v.s. the elastic constant k n . (b) the same data D y are plotted 
v.s. the average coordination number of the pile z. Each plot represents average over around 350 
realizations. 
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oUSkjO (a) Ftictionless ^.U.y) (b) Friction less 




FIG. 7: (Color online) The spatial distribution of the averaged absolute value of the longitudinal 
response in the frictionless piles (a, b), and in the frictional piles (c, d) with k n = 10 3 [mo<?/cro] and 
10 6 [mo<7/<To]- The plots are averaged over a few hundreds realizations. 
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FIG. 8: (Color online) The spatial distribution of the averaged displacement-displacement response 
function in the frictionless piles with k n = 10 3 [mo<?/cro] (a) and 10 6 [mo<7/<7o] (b). The plots are 
averaged over a few hundreds realizations. 
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(a) Friclional 



(b) Friclional 




FIG. 9: (Color online) The spatial distribution of the averaged displacement-displacement response 
function in the frictional piles with k n = 10 3 [mog/<7o] ( a ) an d 10 6 [mo<7/o"o] (b). The plots are 
averaged over a few hundreds realizations. 
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(a) Frictionless 



(b) Frictional 




FIG. 10: (Color online) The contour plots of the longitudinal component of displacement- 
displacement response function A y in the frictionless piles (a) and in the frictional piles (b) with 
the disk elasticity k n = 10 3 , 10 4 , and 10 6 [mog/ao]- 
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